#ifndef __Equation_Solver__
#define __Equation_Solver__

#include<iostream>
#include<limits>
#include<stdexcept>
#include<vector>
#include<algorithm>
using namespace std;

#define MYFUNC(x) double (*(x))(double) 
#define ll long long 
const double epsilon = 10*std::numeric_limits<double>::epsilon();

/*	func类：
	虚函数方法，具体函数通过main定义子类实现；
	原函数默认为identify, 导函数默认数值方法；
	对于高阶导数，默认为1阶，其余需要自行给出；
*/
class func{
public:
	virtual double operator () (double x){
		return x;
	}
	virtual double diff(double x){
		return (operator () (x+epsilon)-operator () (x-epsilon))/(2*epsilon);
	}
	virtual double deriv(double x, int p){
		return (operator () (x+epsilon)-operator () (x-epsilon))/(2*epsilon);
	}
};
/*
	EquationSolver类：
	求解器方法在子类定义，默认为输出常数，意在勉励：
	逸一时，误一世
*/
class EquationSolver{
private:
	func& _target;
	ll times;
public:
	EquationSolver(func& f):_target(f),times(0){}
	double f(double x){return _target(x);}
	double diff(double x){return _target.diff(x);}
	virtual double solve(){	return 114514; }
	void iter(){times++;}
	void print_iters(){std::cout<< "Times of iteration: " << times << std::endl;}
};
/*
	Bis_Solver类:
*/
class Bis_Solver:public EquationSolver{
private:
	double lower;
	double upper;
	double delta;
	ll M; 
public:
	/*
	提供调试类型和使用类型构造函数
	*/
	Bis_Solver(func& f, double lower, double upper, double delta,ll M):EquationSolver(f),lower(lower),upper(upper),delta(delta),M(M){}
	Bis_Solver(func& f, double lower, double upper):EquationSolver(f),lower(lower),upper(upper),delta(epsilon),M(100){}
	double solve(){
		ll iters;
		double c = 0;
		for(iters=1;iters<=M;iters++){
			iter();
			c = (lower + upper)/2;
			if(f(c)==0&&upper-lower<delta) return c;
			if(f(c)*f(lower)>0) lower = c;
			else upper = c;
		}
		return c;
	}
};

class Newton_Solver:public EquationSolver{
private:
	double delta;
	ll M;
	double x_0;
public:
	Newton_Solver(func& f,double x_0,double delta,ll M):EquationSolver(f),delta(delta),x_0(x_0),M(M){}
	Newton_Solver(func& f,double x_0):EquationSolver(f),delta(epsilon),x_0(x_0),M(100){}
	double solve(){
		double curr = x_0-f(x_0)/diff(x_0);
		double prev = x_0;
		for(ll i=1;i<=M;i++){
			iter();
			if(f(curr) == 0&&(abs(curr-prev)<delta)) return curr;
			prev = curr;
			curr = curr - f(curr)/diff(curr);
		}
		return curr;
	}
};

class Scant_Solver:public EquationSolver{
private:
	double delta;
	ll M;
	double x_0,x_1;
public:
	Scant_Solver(func& f,double x_0,double x_1, double delta,ll M):EquationSolver(f),delta(delta),x_0(x_0),x_1(x_1),M(M){}
	Scant_Solver(func& f,double x_0,double x_1):EquationSolver(f),delta(epsilon),x_0(x_0),x_1(x_1),M(100){}
	double solve(){
		double x=x_0;
		double y=x_1;
		for(ll i=1;i<=M;i++){
			iter();
			if((abs(x-y))<delta) return (x+y)/2;
			double temp = x-f(x)*(x-y)/(f(x)-f(y));
			y = x;
			x = temp;
		}
		return (x+y)/2;
	}
};
/*
	DivDiff函数：
	用于计算一列互不相同实数的差商
*/
double DivDiff(func& f, vector<double> _input){
	vector<double> zeros = _input;
	sort(zeros.begin(),zeros.end());
	vector<double> value;
	for(int i=0;i<zeros.size();i++){
		if(i==0) for(int j=0;j<zeros.size();j++) value.push_back(f(zeros[j]));
		else{
			vector<double> temp;
			for(int j=0;j+i<zeros.size();j++){
				if(zeros[j]==zeros[j+i]){
					double ans = f.deriv(zeros[j],i);
					for(int k=1;k<=i;k++) ans/=k;
					temp.push_back(ans);
				}
				else temp.push_back((value[j+1]-value[j]/(zeros[j+i]-zeros[j])));
			}
			value = temp;
		}
	}
	return value.back();
}
/*
	多项式类
*/
class polynomial{
private:
	vector<double> _coeffi;
	void _nonzeros(){
		if(_coeffi.size()<0) return;
		while(_coeffi.back()==0&&_coeffi.size()>=2) _coeffi.pop_back();
	}
public:
	polynomial(){
		_coeffi.push_back(0);
	}
	polynomial(vector<double> _input){
		_coeffi = _input;
		if(deg()<0) _coeffi.push_back(0);
		_nonzeros();
	}
	polynomial(double x_0,double x_1){
		_coeffi.push_back(x_0);
		_coeffi.push_back(x_1);
		_nonzeros();
	}
	~polynomial(){
		_coeffi.clear();
	}
	double operator () (double x){
		double ans = 0;
		for(int i=_coeffi.size()-1;i>=0;i--){
			ans*=x;
			ans+=_coeffi[i];
		}
		return ans;
	}
	int deg(){
		return _coeffi.size()-1;
	}
	double operator [] (int i){
		if(i>deg()) return -1;
		return _coeffi[i];
	}
	polynomial operator + (polynomial& plus){
		vector<double> _plus;
		int i;
		for(i=0;i<=deg()&&i<=plus.deg();i++) _plus.push_back(_coeffi[i]+plus[i]);
		for(;i<=deg();i++) _plus.push_back(_coeffi[i]);
		for(;i<=plus.deg();i++) _plus.push_back(plus[i]);
		return polynomial(_plus);
	}
	void operator += (polynomial& plus){
		vector<double> _plus;
		int i;
		for(i=0;i<=deg()&&i<=plus.deg();i++) _plus.push_back(_coeffi[i]+plus[i]);
		for(;i<=deg();i++) _plus.push_back(_coeffi[i]);
		for(;i<=plus.deg();i++) _plus.push_back(plus[i]);
		_coeffi = _plus;
		_nonzeros();
		return;
	}
	polynomial operator * (polynomial& time){
		vector<double> _time;
		for(int i=0;i<=deg()+time.deg();i++) {
			double _cof = 0;
			for(int j=0;j<=i;j++){
				if(i-j>time.deg()||j>deg()) continue;
				_cof += _coeffi[j]*time[i-j];
			}
			_time.push_back(_cof);
		}
		return polynomial(_time);
	}
	void operator *= (polynomial& time){
		vector<double> _time;
		for(int i=0;i<=deg()+time.deg();i++) {
			double _cof = 0;
			for(int j=0;j<=i;j++){
				if(i-j>time.deg()||j>deg()) continue;
				_cof += _coeffi[j]*time[i-j];
			}
			_time.push_back(_cof);
		}
		_coeffi = _time;
		_nonzeros();
		return;
	}
	void cons_plus(double x){
		_coeffi[0]+=x;
		return;
	}
	double diff(double x){
		double ans = 0;
		for(int i=_coeffi.size()-1;i>=1;i--){
			ans*=x;
			ans+=_coeffi[i]*i;
		}
		return ans;
	} 
	void print(){
		for(int i=0;i<=deg();i++){
			cout << _coeffi[i] << " ";
		}
		cout << endl;
	}
	void latex(){
		for(int i=0;i<=deg();i++){
			if(i==0) cout << _coeffi[i];
			else{
				if(_coeffi[i]>0) cout << "+";
				if(i==1) cout << _coeffi[i] << "x";
				else cout << _coeffi[i] << "x^" << i;
			}
		}
		cout << endl;
	}
};
/*
	Newton插值类:
	必须输入一个vector与函数用以构造多项式
*/
class Newton_inter{
private:
	func& _f;
	vector<double> _divdiff;
	vector<double> _zeros;
	double _DivDiff(func& f, vector<double> _input){
		vector<double> zeros = _input;
		sort(zeros.begin(),zeros.end());
		vector<double> value;
		for(int i=0;i<zeros.size();i++){
			if(i==0) for(int j=0;j<zeros.size();j++) value.push_back(f(zeros[j]));
			else{
				vector<double> temp;
				for(int j=0;j+i<zeros.size();j++){
					if(zeros[j]==zeros[j+i]){
						double ans = f.deriv(zeros[j],i);
						for(int k=1;k<=i;k++) ans/=k;
						temp.push_back(ans);
					}
					else temp.push_back((value[j+1]-value[j])/(zeros[j+i]-zeros[j]));
				}
				value = temp;
			}
			_divdiff.push_back(value[0]);
		}
		return value.back();
	}
public:
	Newton_inter(func& f,vector<double> input):_f(f),_zeros(input){
		if(input.size()==0) return;
		_DivDiff(f,input);
		sort(_zeros.begin(),_zeros.end());
		return;
	} 
	~Newton_inter(){
		_divdiff.clear();
		_zeros.clear();
	}
	double operator () (double x){
		double ans = 0;
		for(int i=_divdiff.size()-1;i>=0;i--){
			ans*=(x-_zeros[i]);
			ans+=_divdiff[i];
		}
		return ans;
	}
	int deg(){
		return _divdiff.size()-1;
	}
	void print(){
		for(int i=0;i<_divdiff.size();i++) cout << _divdiff[i] << " ";
		cout << endl;
	}
	void latex(){
		for(int i=0;i<_divdiff.size();i++){
			cout << _divdiff[i] << "\\pi_{" << i << "}(x)";
			if(i>=0&& i<_divdiff.size()-1){
				if(_divdiff[i+1] >=0 ) cout << "+";
			}
		}
		cout << endl;
		return;
	}
	polynomial stdrize(){
		polynomial ans;
		for(int i=_divdiff.size()-1;i>=0;i--){
			polynomial temp(-_zeros[i],1);
			//cout << _divdiff[i] << " ";
			//temp.print();
			ans*=temp;
			ans.cons_plus(_divdiff[i]);
		}
		return ans;
	}
};

typedef Newton_inter Hermite_inter;
#endif